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Abstract 

Two  deformable  mirrors  with  finite  conjugate  ranges  are  investigated  for 
compensating  amplitude  and  phase  distortions  due  to  laser  propagation  through  turbulent 
atmospheres.  Simulations  are  performed  based  on  Adaptive  Optics  (AO)  for  an  Airborne 
Laser  (ABL)-type  scenario. 

The  Strehl  ratio,  the  number  of  branch  points  in  deformable  mirror  (DM)  controls, 
and  the  number  of  iterations  to  convergence  are  used  as  figures  of  merit  to  evaluate 
performance  of  the  Sequential  Generalized  Projection  Algorithm  (SGPA)  that  generates 
mirror  commands.  The  Strehl  ratio  and  the  number  of  branch  points  are  plotted  versus 
the  log-amplitude  variance  (also  known  as  the  Rytov  parameter),  the  conjugate  range  of 
the  second  deformable  mirror,  and  the  radii  of  the  deformable  mirrors.  Also,  the  number 
of  iterations  is  plotted  versus  the  Rytov  parameter  and  the  conjugate  range  of  the  second 
deformable  mirror.  The  results  are  ensemble  averages  over  32  realizations  of  the 
scintillated  test  fields  for  each  value  of  the  Rytov  parameter  within  the  test  scenario. 

The  Gaussian  beam  shape  that  optimizes  the  Strehl  ratio  is  determined.  The  least 
squares  two-DM  Strehl,  phase-only  Strehl,  least  squares  phase-only  Strehl,  and 
uncompensated  Strehl  are  also  determined  for  comparison.  Finally,  for  the  Strehl  ratio 
versus  Rytov  parameter  analysis  the  Strehl  is  also  calculated  beyond  the  telescope  by 
propagating  the  pre-compensated  laser  wavefront  back  through  the  phase  screens  of  the 
modeled  atmosphere. 

A  conclusion  is  that  an  AO  transmission  system  with  two  DMs  clearly  improves 
theoretical  performance  (compared  to  a  system  with  one  DM)  in  delivering  energy  on 
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target  through  atmospheric  turbulence.  Also,  placing  the  second  DM  at  a  finite  conjugate 
range  minimizes  energy  lost  outside  the  radius  of  the  first  deformable  mirror  and  thus 
maximizes  the  Strehl  ratio.  Finally,  it  is  also  concluded  that  hidden  phase  contained  in 
the  branch  points  is  critical  to  the  performance  of  the  SGPA  algorithm.  It  is  suggested 
that  non-least  squares  methods  and/or  branch  point  number  constraints  could  reduce 
hidden  phase  effects  to  further  improve  performance. 
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MULTI-CONJUGATE  ADAPTIVE  OPTICS  FOR  THE  COMPENSATION 


OF  AMPLITUDE  AND  PHASE  DISTORTIONS 


1.  Introduction 

A  majority  of  Adaptive  Optics  (AO)  research  to  date  has  focused  on  phase-only 
correction  systems,  where  the  primary  goal  is  the  improved  imaging  of  space  objects  at 
observation  angles  near  zenith.  For  laser  communications  or  high-energy  weapon 
systems  such  as  the  Airborne  Laser  (ABL),  however,  AO  must  also  provide  correction 
for  scintillation  in  the  amplitude  of  the  transmitted  field.  This  thesis  explores  the  two- 
deformable  mirror  (DM),  multi-conjugate  configuration  applied  to  full  optic  conjugation 
in  order  to  correct  for  both  phase  and  amplitude  distortions  in  the  transmitted  beam. 

1.1  Problem  Statement 

High  quality  telescope  mirrors  as  large  as  10  meters  are  available,  but  typically 
they  can  provide  angular  resolution  no  better  than  that  of  a  25  cm  telescope  at  optical 
wavelengths  due  to  atmospheric  turbulence.  The  Hubble  Space  Telescope  (HST)  is  one 
solution  to  this  problem  but  is  very  expensive.  Another  solution  is  adaptive  optics  [1]. 

Adaptive  optics  systems  can  be  broadly  divided  into  two  basic  types,  those 
designed  for  transmission  and  those  designed  for  imaging.  Examples  of  transmitting 
systems  are  tactical  communications  and  high-energy  laser  weapons  (e.g.,  ABL). 

Imaging  systems  are  used  mainly  for  astronomy  or  ground-based  satellite  surveillance. 
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Common  to  all  AO  systems,  however,  are  the  performance-hindering  effects  of 
atmospheric  turbulence. 

For  transmission  systems,  turbulence  can  result  in  degradation  of  the  signal  to 
noise  ratio  or  in  a  decrease  in  energy  on  target.  The  goal  addressed  in  this  thesis  is  the 
evaluation  of  a  deformable  mirror  control  algorithm  for  a  Multi-Conjugate  Adaptive 
Optics  (MCAO)  system  designed  to  compensate  for  strong  turbulence  due  to  propagation 
over  long  horizontal  paths,  where  the  second  deformable  mirror  is  conjugate  to  a  finite 
distance  from  the  telescope’s  collecting  aperture. 

1.2  Document  Organization 

Chapter  two  summarizes  various  background  theories  on  which  this  research  is 
based.  Chapter  three  introduces  the  Fourier-transform-based  phase  screen  generation 
method  used  to  model  atmospheric  turbulence.  Layered  models  together  with  a  wave 
propagation  algorithm  are  used  to  create  test  fields  with  varying  degrees  of  scintillation  in 
the  amplitude  profile.  Statistics  of  both  the  individual  phase  screens  and  the  scintillated 
test  fields  are  calculated  and  compared  with  theory.  Chapter  four  develops  the  algorithm 
used  to  obtain  the  deformable  mirror  commands  and  describes  branch  points  and  least 
squares  reconstruction.  Chapter  five  presents  and  discusses  key  results.  Finally,  chapter 
six  summarizes  the  research  effort  and  provides  recommendations  for  further  research. 
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2.  Background 


2.1  Adaptive  Optics 

It  has  long  been  known  that  the  turbulence  in  the  atmosphere  distorts  images  of 
heavenly  bodies  as  seen  from  earth.  Operation  of  a  single  deformable  mirror  AO 
imaging  system  to  compensate  for  these  distortions  is  conceptually  quite  simple.  Light 
from  a  distant  point  source  is  essentially  a  plane  wave  as  it  enters  the  earth’s  atmosphere. 
Propagating  downward  through  the  atmosphere,  a  spatially  coherent  planar  wave 
encounters  pockets  of  air  that  vary  in  temperature  and  therefore  density,  resulting  in 
slight  variations  in  refractive  index.  Different  portions  of  the  wavefront  are  thus  subject 
to  slightly  different  optical  path  lengths  between  the  top  of  the  atmosphere  and  the 
collecting  aperture  of  a  telescope. 

The  distorted  wavefront  arriving  at  the  telescope  is  often  described  as  being 
wrinkled.  These  wrinkles,  or  distortions,  can  be  monitored  in  real  time  by  wavefront 
sensors.  Popular  designs  include  the  Shack-Hartmann  wavefront  sensor  and  the  shearing 
interferometer.  Imaging  the  pupil  plane  of  the  telescope  onto  a  deformable  mirror  that  is 
controlled  on  a  millisecond  timescale  allows  these  “wrinkles”  to  be  removed,  resulting  in 
images  of  near  diffraction-limited  quality.  A  tip-tilt  mirror  is  also  normally  included  in 
the  system  to  correct  for  overall  tilt  in  the  distorted  wavefront,  thus  removing  what 
appears  as  jitter  to  the  human  eye. 

2.1.1  Brief  History 

Like  many  other  ideas  that  were  realized  for  the  first  time  in  the  twentieth  century 
due  to  technological/manufacturing  advances,  Adaptive  Optics  is  not  particularly  new. 
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Horace  Babcock  was  apparently  the  first  to  suggest,  early  in  the  1950s,  a  system  for 
improved  astronomical  imaging.  His  approach  used  an  electrostatically  controlled  thin 
layer  of  oil  to  introduce  corrective  phase  delays.  In  1957,  a  Russian,  Vladimir  P.  Linnik, 
independently  described  the  same  concept  in  a  Soviet  Journal.  Several  decades  later, 
when  the  space  programs  of  Russian  and  the  U.S.  were  in  full  swing,  the  DOD  took  the 
lead  in  advancing  AO  technology.  The  first  fully  operational  adaptive  optics  system  was 
built  and  installed  in  a  surveillance  telescope  at  Haleakala  Observatory  on  Maui,  Hawaii, 
for  the  purpose  of  imaging  Russian  satellites  launched  during  the  cold  war  [1]. 

2. 1. 2  Single  Deformable  Mirror  System 

A  basic  system  diagram  of  a  single  deformable-mirror  AO  system  applied  to 
imaging  is  shown  below  in  Figure  2. 1 .  The  incoming  distorted  wavefront  is  captured  by 
the  primary  mirror,  re-imaged  onto  the  deformable  mirror,  and  then  passed  on  to  a 
wavefront  sensor.  The  Shack-Hartmann  sensor  actually  measures  wavefront  slopes  over 
an  array  of  sub-apertures.  This  slope  data  is  passed  to  a  computer  that  employs  an 
algorithm  (such  as  least  squares)  to  reconstruct  the  phase  from  the  slopes,  which  is  then 
used  to  provide  control  commands  to  the  deformable  mirror.  The  wavefront  sensor  is 
normally  placed  “downstream”  from  the  deformable  mirror  so  that  the  system  can 
perform  as  a  closed-loop  control  system.  For  observations  near  zenith  the  sensed 
wavefront  contains  primarily  phase-only  distortions,  and  thus  near  diffraction-limited 
results  may  be  obtained. 
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Incoming  Beacon  Field 


Shape  Signals 


Figure  2.1.  Block  Diagram  for  a  Single  Deformable  Mirror  Imaging  System 


2.1.3  Two  Deformable  Mirror  System 

The  addition  of  a  second  deformable  mirror  to  a  transmitting  system  as  shown  in 
Figure  2.2  allows  pursuit  of  full  optic  conjugation.  Such  conjugation  is  achieved  by 
taking  advantage  of  the  fact  that  phase  modulation  followed  by  propagation  results  in 
changes  to  the  light  wave  amplitude.  By  employing  a  deconvolution  type  algorithm,  the 
second  deformable  mirror  can  be  driven  to  a  shape  that  results  in  an  amplitude  field  equal 
to  that  of  the  sensed  field  after  propagation  of  a  laser  from  DM2  to  DM1.  A  two-DM  AO 
system  with  the  second  deformable  mirror  conjugate  to  infinity  (i.e.,  the  far-field)  was 
previously  studied  by  Roggemann  and  shown  to  provide  an  increase  in  on-target  on-axis 
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light  field  amplitude  by  a  factor  of  1.4  to  1.5  as  compared  with  a  one-DM  phase-only 
correction  system  [2], 


Receiving  and 


Field  From  Beacon 


Laser 


Figure  2.2.  Block  Diagram  for  a  Two  Deformable  Mirror  Transmission  System 


2. 1. 4  Figures  of  Merit 

One  of  the  primary  metrics  used  to  assess  the  performance  of  adaptive  optics 
systems  is  the  Strehl  ratio,  which  is  defined  as  the  on-axis  intensity  in  the  far-field 
produced  by  a  system  with  no  aberrations  divided  by  the  on-axis  intensity  of  the  aberrant 
system.  This  ratio  is  expressed  in  terms  of  the  incoherent  imaging  point-spread  function 
(psf), 

(2i 

PsLmouA  0,0) 
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where  with  indicates  the  aberrant  system  and  without  indicates  the  diffraction  limited 
system.  The  amplitude  psf  describes  a  coherent  imaging  system  and  is  simply  the 
Fraunhofer  diffraction  pattern  of  the  exit  pupil, 


=  J dxdy  P{x,y)tx^-i^-{ux+vy) 


2 n 

Xzs 


(2.2) 


whereof  is  a  constant  amplitude,  A  is  the  wavelength  of  light,  z(.  is  the  distance  between 

the  exit  pupil  and  image  planes,  P(x,  y)  describes  the  shape  of  the  exit  pupil,  and  (x,  y) 
and  (w,  v)  are  the  coordinates  in  the  exit  pupil  and  image  planes,  respectively.  The 
amplitude  transfer  function  is  simply  the  Fourier  transform  of  Equation  (2.2),  which  turns 
out  to  be  simply  a  scaled  version  of  the  pupil  function.  Because  the  pupil  function  is 
normally  symmetric  the  amplitude  transfer  function  can  be  expressed  as 

H  {fx’fr)  =  •  (2-3) 

The  Optical  Transfer  Function  (OTF)  describes  an  incoherent  imaging  system  and  is  the 
normalized  Fourier  transform  of  the  squared  modulus  of  the  amplitude  psf, 

|  J  dudv\h(u,v)\  exy[-j2n(fxu+fYv)\ 

{/x’/y)=  •  (2.4) 

J  J  du  dv\h(u,v)^ 


With  the  help  of  Rayleigh’s  theorem,  the  OTF  can  also  be  expressed  as  the  normalized 
autocorrelation  function  of  the  amplitude  transfer  function 


)  — " 


j  J  dxdyP 


2,Zjx  kzjy 
x+ — ~,y+ — — 

r\  ?  ✓  r\ 


yix-^djL.y-^h- 


(2.5) 


J  J  dxdy\P(x,y)f 
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For  an  aberrant  system  the  psf  is  calculated  in  the  same  fashion  after  replacing  P(x,y) 
with  a  generalized  pupil  function  given  by 

P(x,y)  =  P{xty)cx$\ikW  (x,y)]  ,  (2.6) 


where  W  (x,  y)  describes  the  remaining  aberrations  present  in  the  system  after  adaptive 
optics  compensation.  The  OTF  describing  an  aberrant  system  is  thus 


P_  {XzJx ,  XzJY  )★?*  (A  z,fx ,  XzJY ) 

p(0,0)*p*(0,0) 


(2.7) 


The  psf  describing  an  incoherent  imaging  system  can  now  be  expressed  as  the  inverse 
Fourier  transform  of  the  OTF  given  by  Equation  (2.7), 

Psf  (x,y)  =  ?''{&  {fx,fY)}  (2.8) 

where  7'  represents  the  two-dimensional  inverse  Fourier  transform. 

Equation  (2.1)  can  be  directly  applied  to  single  deformable  mirror  AO  imaging 
systems,  but  it  must  be  altered  somewhat  for  application  to  the  2-DM  laser  transmission 
configuration.  In  order  to  reverse  the  effects  of  turbulence,  the  outgoing  laser  wavefront 
must  have  an  amplitude  profile  identical  to  that  of  the  sensed  beacon  field  and  a  phase 
profile  equal  to  the  conjugate  of  the  sensed  wavefront.  The  far-field  on-axis  intensity  is 
clearly  proportional  to  the  degree  to  which  phase  conjugation  is  successfully 
accomplished,  or 

Hl-.WMf  .  (2-9) 
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where  U0  (r)  =  \  (r)e“t>"^  describes  the  field  of  the  outgoing  laser  in  the  pupil  plane  and 

Ub  (r)  =  Ab  (r)eh^  describes  the  sensed  beacon  field  in  the  pupil  plane.  In  the  best-case 
scenario  the  field  of  the  outgoing  laser  is 

Uo,  ideal  (r)  =  kAb  (r)e'h{r]  ,  (2. 1 0) 

where  A:  is  a  constant.  The  maximum  far-field  on-axis  intensity  is  now 

I  =\[drkAb(r)Ab(r)\  =  f  dr  k2Ab2  (r)j  d/ Ab2  (/)  .  (2.11) 

max 

The  normalized  far-field  on-axis  intensity  is  then 

\$drU0(r)U,(rf 

The  first  factor  in  the  denominator  represents  the  energy  contained  in  the  field 
transmitted  by  the  laser  after  modulation  by  DM2,  propagation  from  DM2  to  DM1,  and 
modulation  by  DM1.  Since  the  ffee-space  propagation  represented  by  Tz  [•]  conserves 

energy,  the  first  factor  in  the  denominator  equals  the  energy  in  the  original  laser  field 
(assuming  no  loss  from  the  deformable  mirrors),  i.e., 

Energy laser  =\drU,(r) U*  ( r)  =  \dr  A,2  (r)  .  (2. 1 3) 


Equation  (2.12)  is  finally  rewritten  as 


2.2  Propagation  of  Optical  Radiation 


2.2.1  Angular  Spectrum  and  the  Propagation  Transfer  Function 

The  propagation  of  optical  fields  through  homogeneous  media  can  be  formulated 
to  allow  the  application  of  linear  systems  theory,  which  is  an  invaluable  tool  for  analysis. 
A  wave  traveling  in  the  positive  z  direction  is  described  in  any  plane  perpendicular  to  its 
path  by  its  angular  spectrum,  which  is  simply  the  Fourier  transform  of  the  field  with 
direction  cosine  terms  substituted  for  the  spatial  frequency  variables.  This  formulation 
indicates  that  any  field  can  be  described  in  terms  of  a  weighted  sum  of  plane  waves 
traveling  in  various  directions.  If  a  field  is  known  at  z  =  0,  then  the  effects  of  propagation 
through  the  homogeneous  media  are  described  by  the  transfer  function  of  wave 
propagation 


*(/,./,)  = 


exp  j2n  j^!-(A/J2-(A4)2  Jfx 


+f  <■ 


(2.15) 


otherwise 


such  that  the  field  after  traveling  a  distance  z  is 

A(fx,fy;z)  =  A(fx,fy;0)H(fx,fy)  ,  (2.16) 

where  A^fx,fy,  0)  is  the  two-dimensional  Fourier  transform 


A(fxJy‘,0)=  J  ^(hdyU(x,y,0)e^\-iln{fxx+fyy)\  ,  (2.17) 

and/r  and  fY  are  the  spatial  frequency  variables  in  the  jc  and  y  directions,  respectively. 
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2.2.2  Fresnel  Approximation  to  the  Angular  Spectrum  Propagator 

A  simple  method  for  deriving  the  Fresnel  approximation  to  the  propagation 
transfer  function  applies  a  binomial  expansion  and  approximation  to  the  exponent  of 

Equation  (2.15).  The  binomial  expansion  of  >/l  -  x  is 


Using  the  first  two  terms  of  the  expansion,  the  exponent  is  simplified,  i.e., 

VMV,)2-(vJ  “i • 

The  Fresnel  approximation  to  the  angular  spectrum  propagator  is  thus 

H(f  /)={ exp,fe exp \rinXz U2 + fy )]  I1/, I < 1  and  \Xf* | <  1 

0  otherwise 


(2.18) 


(2.19) 


(2.20) 


2.3  Statistics  and  Random  Processes 

For  linear  time-invariant  (LTI)  systems  and  deterministic  signals,  the  input,  x(t), 
and  output,  y(t),  are  related  by 

y(t)  =  x(t)®h(t)  ,  (2.21) 

where  h(t)  is  the  impulse  response  of  the  system  and  ®  is  the  convolution  operator.  The 
input-output  relation  is  expressed  by  straightforward  multiplication  in  the  frequency  (i.e., 
Fourier  Transform)  domain 

Y(f)  =  X{f)H(f)  ,  (2.22) 
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where  Y(f),  X(f),  and  H(f)  are  the  Fourier  transforms  of  y(t),  x(t),  and  h(t),  respectively. 
Under  certain  conditions  a  similar  relationship  holds  for  stochastic  signals.  For  Wide 
Sense  Stationary  (WSS)  random  processes,  the  power  spectral  density  of  the  input  and 
output  of  a  LTI  system  are  related  by 

’  <2-23) 

where  Sx(f)  and  Sy(f)  are  the  power  spectra  of  the  input  and  output,  respectively,  and  H(f) 
is  again  the  Fourier  transform  of  the  system  impulse  response.  Another  special 
relationship  exists  between  the  autocorrelation  function  and  power  spectral  density  of  a 
WSS  random  process: 

S(/)  =  ?{r(l)}  .  (2.24) 

In  general,  the  spatial  autocorrelation  function  of  a  real-valued  process  <p  (r)  is 

r*(ri  =  .  (2.25) 

Several  assumptions  greatly  simplify  this  expression.  First,  if  the  process  is  assumed  to 
be  homogeneous,  then  the  statistics  describing  it  are  independent  of  location,  and  the 

autocorrelation  is  only  a  function  of  the  separation  vector  r  =  rx-r2,  i.e., 

{r)  =  E{<l>(ri)<t,{n-r)}  •  (2.26) 

Second,  if  the  process  is  also  isotropic,  then  the  statistics  (i.e.,  properties  of  the  medium) 
do  not  depend  on  direction,  and  the  autocorrelation  is  a  function  only  of  the  magnitude  of 
the  separation  vector. 

Subsequent  to  propagation  through  atmospheric  turbulence,  the  variance  of  an 
uncompensated  distorted  wavefront  observed  in  a  circular  aperture  is 
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cr  =1.0299 

uncomp  * 


fD^ 


Vr*J 


(2.27) 


where  D  is  the  aperture  diameter  and  r0  is  the  Fried  parameter,  which  is  related  to 
turbulence  strength  [3].  Subsequent  to  piston  and  tilt  removal,  the  variance  is 


<ncomp  =  0.134 


rD^ 


(2.28) 


The  corresponding  expression  for  a  square  aperture  is 


=0-1748 


vr°  J 


(2.29) 


where  L  is  the  length  of  a  side  of  the  aperture. 

Another  statistical  quantity  used  frequently  to  describe  turbulence  is  the  structure 
function 


^(r1,r2)  =  ^(r1)-0(r2)]2}  .  (2.30) 

For  homogeneous  media,  the  structure  function  is 

•  (2-31) 

By  expanding  the  argument  of  the  expectation  operator,  Equation  (2.31)  can  be  rewritten 
as 


D*  {r  )  =  E  {  *  (ri  /  -  (*i )  0  fa  ~  r )  +  0  (*i  -  rf  J} 


The  homogeneity  assumption  (i.e.,  E[<j)2  {rx )]  =  E\j)2  (r2)])  allows  further 
simplification: 


(2.32) 
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(2.33) 


Finally,  using  the  definition  of  the  autocorrelation  function  given  by  Equation  (2.26),  the 
structure  function  is 

fl,(?)  =  2[r#(°)-r,(r)]  .  (2.34) 

For  plane  waves  the  phase  structure  function  for  Kolmolgorov  turbulence  is  [4:34] 

(  r  f 3 

Df(r)  =  6.88  —  ,  (2.35) 

v°  J 

where  r  is  the  magnitude  of  the  radial  position  vector  in  a  plane  perpendicular  to  the 
direction  of  propagation. 

2.4  Atmospheric  Turbulence  Modeling 

One  of  the  effects  of  turbulence  on  the  propagation  of  light  through  the 
atmosphere  is  image  blurring,  which  is  due  primarily  to  the  negative  effects  of  turbulence 
on  phase.  For  a  plane  wave  entering  the  Earth’s  atmosphere  with  an  arbitrarily  large 
spatial  coherence,  propagation  over  a  distance  L  through  the  turbulent  atmosphere  results 
in  degraded  spatial  coherence  in  the  aperture  plane  of  a  collecting  telescope.  This 
degradation  is  measured  by  the  transverse  coherence  length  [5], 

[O  Q 1  "1-3/5 

~Yk2\ldz  cn  00  .  <  P  <  L0  (2.36) 

3  44  (of3  cl  V'5 

Po=  —k2  j  \QdzCl{z)  >  P< 4  »  (2.37) 

.  V°) 
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where  p  is  the  magnitude  of  the  radial  position  vector  in  the  plane  of  the  aperture,  po  is 
the  separation  distance  beyond  which  two  points  on  the  wavefront  in  the  aperture  are 
uncorrelated,  lo  is  the  inner  scale  which  corresponds  to  the  smallest  eddy  size,  L0  is  the 
outer  scale  which  corresponds  to  the  largest  eddy  size,  and  k  is  the  optical  wavenumber 
(InlX).  Two  key  relationships  should  be  noted  based  on  Equations  (2.36)  and  (2.37). 
First,  the  transverse  coherence  length  is  proportional  to  X6/s,  and  therefore  image 
degradation  due  to  turbulence  is  somewhat  less  in  the  infrared  region  than  in  the  visible. 
Second,  po  is  proportional  to  the  integrated  turbulence  strength.  The  factor  C2  (z)  is  the 

refractive  index  structure  constant  and  is  a  measure  of  the  strength  of  turbulence,  but  it  is 
not  truly  constant  and  varies  with  time,  geographic  location  on  the  earth,  and  altitude. 
Since  much  of  the  turbulence  occurs  at  lower  altitudes,  major  observatories  around  the 
world  are  built  on  high  mountaintops  to  effectively  reduce  the  integrated  turbulence 
strength. 

A  much  more  common  parameter  that  describes  the  negative  effects  of  turbulence 
on  spatial  coherence  is  the  atmospheric  coherence  length  (also  known  as  the  Fried 
parameter)  r0,  which  is  a  translation  of  the  transverse  coherence  length  from  the  aperture 
plane  to  the  focal  plane.  This  parameter  was  defined  in  1966  by  David  Fried  to  be  the 
largest  aperture  within  which  the  total  root  mean  square  (rms)  wavefront  irregularity  is 
less  than  one  radian  (or  7J2ti)  [1],  The  Modulation  Transfer  Function  (MTF)  in  the  focal 
plane  [6]  is 

MTF(v)  =  exV[-3M(Xfv/r0f]  ,  (2.38) 
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where /is  the  focal  length  of  the  system,  v  is  the  spatial  frequency,  and  r0  is  the 
atmospheric  coherence  length,  which  is  related  to  the  transverse  coherence  length  by  r0  = 
2.1  po.  The  atmospheric  coherence  length  is  thus 


r0  =1.67691 


.£  -1-3/5  _ 

k2)^  dzC2(z)  ,yfXL<p<L0 


(2.39) 


r0  =1.5167 


-i-3/5 


L 


v°y 


tdzCl(z) 


,  P<h 


(2.40) 


The  effective  MTF  of  an  optical  system  looking  through  a  turbulent  atmosphere  is  simply 
the  product  of  the  system  MTF  and  the  MTF  describing  the  atmosphere.  Therefore 
turbulence  has  the  effect  of  a  low  pass  filter.  The  elimination  of  higher  frequencies  then 
results  in  blurring  of  the  image. 

The  amplitude  and  intensity  of  the  scintillated  field  are  usually  modeled  as  log¬ 
normal  random  variables.  The  log  amplitude  covariance  function  [5:181]  is 

K2  (L-7)) 


Cx(p,L)  =  4 n2k2 £  drj £  dK kJ0  (rep) sin2 


2k 


<I>(K,ri)  ,  (2.41) 


where  p  is  the  magnitude  of  the  radial  position  vector,  L  is  the  propagation  length,  k  is  the 
optical  wavenumber  (2n/X),  k*  is  the  three-dimensional  spatial  wavenumber,  T]  is  the 
variable  of  integration  along  the  propagation  path,  and  <P(K,ri)  is  the  power  spectral 
density  (PSD)  of  turbulence.  The  simplest  PSD  used  to  describe  fluctuations  in  the  index 
of  refraction  in  the  atmosphere  is  the  Kolmogorov  spectrum  [4:30], 

<D„0c)  =  0.033C„VI1/3  ,  (2.42) 
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where  C2  is  the  refractive  index  structure  constant  and  Ki s  the  three-dimensional  spatial 

wavenumber.  The  spectrum  of  phase  fluctuations  is  expressed  in  terms  of  r0  and  spatial 
frequency  [4:35]: 

%  {fxJy)  =  0.023r0-5/3  (f2  +  f2  )"‘V6  .  (2.43) 

Several  other  forms  of  the  PSD  of  turbulence  have  been  introduced  [5:174], 
mostly  for  reasons  of  mathematical  convenience.  The  following  PSD  [7]  allows  for 
inclusion  of  a  finite  inner  scale: 


<£„(x:)  =  0.033C„V11/3  exp 


KZ 


(2.44) 


where  K2m  =5.92 /l0  and  lo  is  the  inner  scale.  The  von  Karman  spectrum  [5:174]  allows 
for  inclusion  of  a  finite  outer  scale, 


On(KT)  =  0.033C2  (k2  +  K2)-11'6  , 


(2.45) 


where  K0  =  2n/L0  and  Lo  is  the  outer  scale.  Finally,  the  modified  von  Karman  spectrum 
[5:174]  allows  for  both  finite  inner  and  outer  scales, 


(k)  =  0.033C2  (k2  +  k2  y' 1/6  exp 


v  K"  ) 


(2.46) 


An  expression  for  the  modified  von  Karman  spectra  that  describes  phase  fluctuations  is 


%{fx,fr)  =  0.023rc 


-5/3 


f2+f!  + 


(  1  \2  'l 


\  L<>  J 


11/6 


exp 


~(f2X+fy) 


\2\ 


0.942 


.(2.47) 


Returning  to  the  discussion  of  scintillation  statistics,  the  log-amplitude  variance  is 
expressed  in  terms  of  the  propagation  length  by  assuming  the  Komologorov  spectrum  for 
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the  PSD  and  evaluating  the  covariance  function  (Equation  (2.41))  at  p  =  0.  The  required 
integration  is  performed  numerically  with  the  following  variable  substitutions: 


a  = , 


i-n 
'  2k 


K 


da-. 


1-7] 


2k 


dK 


(2.48) 


which  results  in  the  simplified  covariance  function 

o2  ( L )  =  (0.132)tt2  (2)~5/s  kVG^dr\  (L-rjf6  C2 (rj) daa~$l$  sin2  (a2)  .  (2.49) 


The  integral  over  a  is  now  performed  numerically,  and  to  four  decimal  places  it  is  equal 
to  0.7701 .  Combining  constants  in  front  of  the  remaining  integral  results  in  a  final 
expression  for  log-amplitude  variance: 

<%!(£)  =  (0.5631 )k^^dn(L-tlfC‘M  ■  (2.50) 

Equation  (2.50)  describes  the  variance  of  the  log  of  the  amplitude  in  a  plane  wave  after 
propagating  a  distance  L.  For  the  log-amplitude  variance  of  an  optical  field  at  the  earth’s 
surface  due,  for  example,  to  an  observed  star,  the  point  7]  =  0  corresponds  to  the  location 
at  which  the  plane  wave  enters  the  earth’s  atmosphere  and  the  point  7]  =  L  is  determined 
by  the  thickness  of  the  atmosphere  traversed.  This  equation  can  also  be  expressed  so  that 
the  zero-point  corresponds  to  the  location  of  a  ground-based  observer.  In  this  case  the 
integration  is  carried  out  from  the  ground  up  to  a  desired  altitude: 

<T=(I)  =  (0.5631)i7',|‘^(^f,C'(/J)  .  (2.31) 

The  important  physical  interpretation  of  the  expressions  for  log-amplitude  variance  is  that 
for  a  given  point  along  the  integration  path  through  the  atmosphere,  the  strength  of 
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turbulence  is  weighted  by  a  factor  proportional  to  the  remaining  distance  to  the 
observation  point  on  the  ground.  For  astronomical  imaging  at  near  zenith,  mild 
turbulence  is  often  modeled  by  a  single  phase  screen  at  the  telescope  pupil.  For  extended 
turbulence  that  results  in  scintillation  in  addition  to  a  wrinkled  phase  profile,  the  phase  is 
distorted  at  some  distance  from  the  pupil  (or  target),  which  means  that  the  individual 
tilted  portions  of  wavefront  eventually  interfere  with  each  other.  The  resulting 
destructive  and  constructive  interference  is  termed  scintillation. 

2.5  Sampling  Theory 

There  are  two  critical  relationships  that  involve  sampling  and  the  Discrete  Fourier 
Transform  (DFT).  First,  for  N  samples  of  a  signal  in  the  spatial  (or  temporal)  domain 
separated  by  Ax,  the  frequency  domain  signal  obtained  by  taking  the  DFT  is  also  repeated 
every  N  samples,  or  1/Ax  in  terms  of  frequency.  Conversely,  N  samples  in  the  frequency 
domain  separated  by  A£,  correspond  to  a  spatial  domain  signal  which  repeats  every  N 
samples,  or  1/A£  in  terms  of  space  (or  time)  units.  These  relationships  are  illustrated  in 
Figure  2.1. 
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Spatial  Domain 


Spatial  Frequency  Domain 


Figure  2.3.  Sampling  and  the  Discrete  Fourier  Transform 


The  other  critical  relationship  pair  relates  the  spacing  in  one  domain  to  the  sampling 
parameters  in  the  other: 


which  is  also 


= 


1 

NAx  ’ 


Ax  = 


1 

NA% 


(2.52) 


(2.53) 


The  matrices  used  in  the  propagation  algorithms  are  setup  as  in  Figure  2.4.  Here 
values  inside  the  matrix  indicate  pixel  distance  from  the  origin.  The  Fresnel 
approximation  to  the  free  space  propagation  transfer  function  (Equation  (2.20))  is  used 
for  all  simulations.  Note  that  to  avoid  aliasing  of  the  propagator,  for  a  generic  phasor  the 
separation  between  samples  of  the  phase  angle  must  be  less  than  or  equal  to  n.  The  phase 
angle  of  the  Fresnel  propagator  is  a  quadratic  function  of  spatial  frequency,  and  thus  the 
largest  difference  between  adjacent  pixels  occurs  at  the  four  comers  of  the  matrix. 
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N±  2 
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N 


Figure  2.4.  Propagation  Matrix  Setup 


Taking  the  lower  right  comer  as  an  example,  the  difference  between  pixel  (N,  N)  and 
pixel  (N-l,  N-l)  must  satisfy  the  previously  stated  requirement.  This  outermost  diagonal 
distance  clearly  contains  the  greatest  adjacent  change.  The  squared  spatial  frequency  in 
any  one  of  the  four  outermost  comers  is,  in  general, 


(2.54) 


The  squared  spatial  frequency  in  the  immediately  adjacent  pixel  along  the  diagonal  is 


(N-2 


-1  + 


(N-2 


Lv 


')'] 


(2.55) 
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Taking  the  difference  between  Equations  (2.54)  and  (2.55),  the  maximum  difference  in 
squared  spatial  frequency  is  now 


AK-2=2(A£) 


(2.56) 


The  condition  imposed  on  the  complex  exponential  describing  Fresnel  propagation  is 


A0=^z2(A|)2 


\(N-2i 

f  JV  -  4  Y 1 

-  + 

LI  2  J 

1  2  JJ 

<K 


(2.57) 


Using  equation  (2.52),  the  aliasing  constraint  is 


2Az(7V-3)<i 

(NAxf 


(2.58) 


If  N  is  much  larger  than  3  then  the  (N-3)  term  in  the  numerator  can  be  replaced  by  N, 
which  allows  Equation  (2.58)  to  be  expressed  in  a  simpler  form: 

2  Xz 


N>- 


(Ax)2 


(2.59) 


Here  N  is  normally  at  least  128  for  propagation  simulations  such  as  those  reported  here, 
and  thus,  the  above  assumption  is  valid. 

Note  that  Equation  (2.59)  can  be  interpreted  as  a  requirement  that  N  be  greater 
than  the  number  of  pixels  that  would  fit  into  the  main  lobe  of  the  Fraunhofer  diffraction 
pattern  resulting  from  a  square  aperture  of  dimensions  equal  to  the  pixel  spacing  Ax. 
Along  one  dimension  in  the  plane  containing  such  a  diffraction  image,  the  width  of  the 
central  intensity  lobe  is 

XR 

x  =  —  ,  (2.60) 

a 
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where  the  x  and  y  axes  are  in  the  plane  of  the  diffraction  pattern  located  a  distance  z  from 
the  aperture,  R  is  the  straight-line  distance  from  the  center  of  the  aperture  (atz  =  0)  to  the 
point  in  the  diffraction  pattern,  and  a  is  the  dimension  of  the  square  aperture.  For  small 
angles  between  the  vector  defined  by  R  and  the  optic  axis,  Equation  (2.60)  clearly 
simplifies  to 


(2.61) 


which  is  easily  recognized  as  the  distance  from  the  origin  to  the  first  intensity  null  in  the 
diffraction  pattern  from  a  narrow  slit  of  width  a.  Comparing  Equations  (2.61)  and  (2.59), 
the  aliasing  constraint  is  interpreted  as  a  requirement  that  the  number  of  pixels  along  one 
dimension  of  the  propagator  matrix  must  be  at  least  equal  to  the  number  of  equally  sized 
pixels  that  would  fit  into  the  diffraction  pattern  created  by  a  single  pixel  located  at  a 
distance  z  from  the  propagator. 
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3.  Phase  Screen  Generation 
3.1  Creation  of  Test  Fields 

The  atmosphere  is  often  modeled  by  a  set  of  discrete  turbulence  layers  [8:66-67]. 
This  layered  approach  is  used  to  create  scintillated  test  fields,  which  are  then  used  to 
evaluate  the  performance  of  the  two  deformable  mirror  phase  conjugating  system.  The 
log-amplitude  variance  is  employed  to  drive  the  other  parameters  of  the  layered  model. 

A  discrete  version  of  Equation  (2.51)  is  needed  to  create  the  model: 

ol{L)  =  0.563\k1/62^6C2„iAzi  ,  (3.1) 

i=i 

where  Mis  the  number  of  layers  in  the  model,  L  is  the  thickness  of  the  atmosphere  being 
traversed,  z,  and  A zi  are  the  height  and  thickness  of  a  specific  layer,  respectively,  and  C2 

describes  the  strength  of  turbulence  in  a  given  layer.  To  further  simplify  this  expression 
the  strength  of  turbulence  in  each  layer  is  expressed  as  a  fraction  of  the  total  integrated 
turbulence, 

WiIc=C2n,Azi  ,  (3.2) 

where  Ic  is  the  total  integrated  value  of  C2  and  Wj  is  a  weighting  factor  corresponding  to 

M 

the  relative  strength  of  turbulence  in  a  given  layer  such  that  =1 .  Equation  (3.2) 

i=\ 

can  now  be  employed  to  rewrite  Equation  (3.1)  as  an  expression  for  Ic  as  a  function  of 
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M<)= 


M 

(0.5631  )k1/6^jWizf 

i= 1 


(3.3) 


The  same  approach  may  be  used  to  modify  Equation  (2.39): 

r% 


r0  =  1.6769 


M 

ei^w, 


i=l 


,  JXL<p<L0 


(3.4) 


Using  the  fact  that  the  weight  vector  W,  sums  to  one  and  substituting  Equation  (3.3)  for 
the  total  integrated  turbulence,  the  atmospheric  coherence  diameter  may  be  expressed  as 
a  function  of  log-amplitude  variance: 

-|3/5 


K)  = 


1.1881 


a Jk 


1  M 

1=1 


,  y[XL<p<L0  . 


(3.5) 


For  the  turbulence  strength  profile  given  below,  values  of  r0  and  C2  are  obtained 
based  on  <72  values  of  0.1,  0.2  ...  1.0.  Once  these  values  are  produced  a  series  of 
scintillated  test  fields  are  created  using  an  outer  scale  value  of  L0  =  «> .  The  phase 

) 

screens  through  which  a  wavefront  is  propagated  to  created  scintillated  fields  are  created 
using  the  well-known  FFT  approach  [9]  and  the  Kolmolgorov  spectrum.  The 
Kolmolgorov  PSD  is  obtained  from  the  von  Karman  PSD  given  by  Equation  (2.47)  by 
letting  the  outer  scale  (L0)  go  to  infinity  and  the  inner  scale  (/«)  go  to  zero, 

d>/^77)  =  0.023r0-5/3(f +772)'n/6  ,  (3.6) 

where  £  and  t]  are  the  spatial  frequency  variables  along  x  and  y,  respectively.  The 
atmospheric  phase  screens  are  generated  by  filtering  complex  white  Gaussian  noise  with 
the  square  root  of  the  Kolmorgorov  power  spectrum,  which  serves  to  force  the  power 
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spectrum  (and  autocorrelation)  to  match  that  of  the  Komolgorov  spectrum.  The  phase 
screens  are  generated  numerically  with  the  aid  of  the  DFT, 

0  {nAx,  mAx)  =  Re{.DFT  { JP{p^q)  [a  (p,  q)  +  ib  ( p ,  tf)]}} 

or  (3.7) 

0  ( nAx,  mAx)  =  Im  j  DFT  ^Jp[p^q)  [a  ( p,q)+ ib  (p,  q )]}} 


where  Ax  is  the  phase  screen  sample  spacing,  P (p,q)  is  the  power  associated  with  the 
spatial  frequency  £  -pAt,  and  rj  =  qAt,,  is  the  sample  spacing  in  Fourier  (frequency) 
space  (see  Equation  (2.52)),  and  a  and  b  are  pseudo  random  numbers  generated  from  a 
Gaussian  distribution  with  zero-mean  and  unit- variance.  The  power  associated  with  each 
sample  in  Fourier  space  is 

P(p,q)  =  <S>t(pA%,mAZ)A%2  (3.8) 


Using  Equations  (3.6),  (3.8),  and  (2.52),  the  phase  screen  is  now  described  by 


(  v6 

0(«Ax,wAx)  =  0.1517  ^  Re|l)FT|(^2+772)  11/12 [a(p,9)  +  j7>(^,^)]|| 


/  y6 

(f)(nAx,mAx)  =  0.1517  y  I  Imjz>FY j(£2  +T]2)  Uln[a{p,q)  +  ib(p,q)^  . 


3.1. 1  NOP  Test  Scenario 

In  support  of  the  ABL  program,  test  data  have  been  collected  recently  at  North 
Oscura  Peak  (NOP)  in  New  Mexico.  A  “beacon”  laser  was  fired  from  an  aircraft  to  a  site 
atop  NOP.  Another  laser,  modified  using  adaptive  optics  driven  by  the  sensed  beacon 
field,  then  was  fired  back  toward  the  aircraft.  The  data  collected  from  the  NOP 
experiments  may  represent  the  actual  engagement  conditions  anticipated  for  ABL.  Three 
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different  aircraft-NOP  radial  distances  were  used  for  collecting  data:  35,  42,  and  54  km. 

For  the  purpose  of  the  NOP  simulations  reported  here,  a  propagation  distance  of  50  km  is 

used.  The  turbulence  experienced  along  that  path  is  modeled  using  ten  phase  screens. 

For  the  simulation  scenario,  the  screen  furthest  from  NOP  is  placed  at  a  horizontal 

distance  of  50  km,  and  the  remaining  screens  are  spaced  evenly  between  the  peak  (0  km) 

and  the  final  screen.  The  strength  of  turbulence  is  assumed  to  be  equal  at  each  of  the 

screens,  resulting  in  the  following  weight  vector: 

W  =[~—  -L  -L  j-  -L  -L  J_  J_  _l  j_1  n  im 

"NOP  LlO  10  10  10  10  10  10  10  10  10  J  • 

Also,  an  optical  wavelength  of  1  pm  is  used  for  all  simulations,  which  is  representative  of 

the  ABL  system.  Table  1  summarizes  the  NOP  test  scenario  simulation  parameters. 


Table  1.  NOP  Test  Scenario  Simulation  Parameters 


C2  (xl0~12  m“2/3) 

r0  (cm) 

2Afr/  Ax2  (pixels) 

0.1 

0.42488 

30.884 

16.77 

0.2 

0.84977 

20.376 

38.54 

0.3 

1.27465 

15.976 

62.69 

0.4 

1.69953 

13.443 

88.54 

0.5 

2.12442 

11.758 

115.72 

0.6 

2.54930 

10.540 

144.03 

0.7 

2.97418 

9.609 

173.29 

0.8 

3.39907 

8.869 

203.41 

0.9 

3.82395 

8.264 

234.29 

1.0 

4.24883 

7.758 

265.86 
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3.2  Comparison  of  Simulated  Phase  Screen  Statistics  to  Theory 

Several  statistics  of  the  FFT-generated  random  phase  screens  are  measured  and 
compared  to  theory  in  the  following  section.  For  all  of  the  following  statistical  analyses, 
a  value  of  10  cm  is  used  for  the  Fried  parameter,  and  the  sampling  interval  is  set  to  ~rQ 

(2.5  cm).  A  sample  screen  is  plotted  as  a  three-dimensional  surface  in  Figure  3.1  and  as  a 
bitmap  in  Figure  3.2. 


Pixels 


128 


Pixels 


Figure  3.1.  Surface  Plot  of  One  Realization  of  a  Random  Phase  Screen  (Using 
Kolmolgorov  Spectrum  with  Ax  =  10  cm  and  r0  =  2.5  cm) 
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Pixels 


Figure  3.2.  Bit-map  of  One  Realization  of  a  Random  Phase  Screen  (Using 
Kolmolgorov  Spectrum  with  Ax  =  10  cm  and  r0=  2.5  cm) 


3.2.1  Phase  Variance  (Tilt  Removed)  vs.  Aperture  Size 

In  order  to  verify  correct  operation  of  the  FFT-based  phase  screen  generator, 
several  realizations  were  created  to  measure  statistical  properties.  The  mean  phase 
variance  versus  square  aperture  width  (normalized  to  r0)  is  plotted  in  Figure  3.3  for  100 
phase  screen  realizations  generated  using  the  Kolmolgorov  spectrum.  The  mean  phase 
variance  versus  circular  aperture  diameter  (normalized  to  r0)  is  also  plotted  in  Figure  3.4 
for  100  phase  screen  realizations  generated  using  the  Kolmolgorov  spectrum.  In  both 
cases  the  variance  is  calculated  for  a  set  of  increasing  aperture  dimensions  by  extracting 
the  pixels  captured  by  a  given  aperture  (mask),  removing  piston  and  Zemike  tilt,  and 
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calculating  the  variance  of  their  values.  One  hundred  vectors  containing  variance  versus 
aperture  size  are  thus  created  and  averaged. 


Figure  3.3.  Phase  Variance  vs.  L/ro  for  a  Square  Aperture 


Figure  3.4.  Phase  Variance  vs.  d/ro  for  a  Circular  Aperture 
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3.2.2  Autocorrelation 

The  mean  autocorrelation  of  the  phase  is  plotted  in  Figure  3.5  for  100  phase 
screen  realizations  generated  using  the  Kolmolgorov  spectrum.  The  autocorrelation  of  a 
single  realization  tends  to  deviate  significantly  from  the  smooth  surface  shown  in  Figure 
3.5  near  the  outer  portions  of  the  discretely  sampled  surface  due  to  the  fact  that  fewer 
data  points  are  available  to  average  as  the  magnitude  of  the  separation  vector  r 
approaches  one-half  the  diagonal  distance  across  the  matrix. 


Figure  3.5.  Phase  Autocorrelation  Function 


3.2.3  Phase  Structure  Function  vs.  d/ro 

The  mean  phase  structure  function  is  plotted  in  Figure  3.6  as  a  three-dimensional 
surface  (using  Equation  (2.34))  for  100  phase  screen  realizations  generated  using  the 
Kolmolgorov  spectrum.  The  averaged  cross  section  of  this  surface  is  displayed  in  Figure 
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3.7.  One  of  the  documented  reasons  for  the  disparity  between  the  theoretical  and  actual 
phase  structure  functions  is  the  loss  of  low  frequency  content  due  to  sampling  of  the 
Kolmolgorov  spectrum  [10].  The  lowest  frequency  in  the  sampled  spectrum  is  l/x,  where 
x  the  physical  width  represented  by  an  array  of  samples.  This  negative  effect  can  be 
offset  to  some  degree  by  increasing  the  number  of  samples  (for  a  constant  spacing),  but 
this  remedy  obviously  imposes  greater  computational  requirements. 


Figure  3.6.  Mean  Phase  Structure  Function  (Surface  Plot) 
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Figure  3.7.  Mean  Cross-Section  of  the  Phase  Structure  Function 

3.3  Comparison  of  Simulated  Scintillated  Amplitude  Field  Statistics  to  Theory 
3.3.1  Log-Amplitude  Variance 

The  mean  log-amplitude  variance  (Rytov  parameter)  values  for  an  ensemble  of  32 
independent  realizations  are  displayed  in  Figure  3.8.  The  saturation  of  the  actual  Rytov 
parameter  for  the  simulated  fields  agrees  with  empirical  data  [5: 1 86-7]. 
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4.  Sequential  Generalized  Projection  Algorithm  (SGPA) 

4. 1  Physics-Based  Algorithm  Development 

Following  the  earlier  discussion  on  wavefront  propagation,  the  operator  Tz  [•] 
indicates  Fresnel  propagation  (given  by  Equation  (2.20))  of  a  field  over  a  distance  z,  and 
T*  [•]  indicates  propagation  in  the  opposite  direction.  The  new  field  at  the  plane  defined 
by  z  is  then 

C/,(7;')  =  r,[i7(,(F)]  =  ?-‘{?[l/„(F)]eM-'’}  ,  (4.1) 

where  7  denotes  the  two-dimensional  Fourier  transform,  z  is  the  distance  of  propagation 
between  the  two  deformable  mirrors,  T  is  the  position  vector  in  the  plane  at  z  =  0,  and  r' 
is  the  position  vector  in  the  plane  at  z  =  z.  The  on-axis  far-field  intensity  of  the  laser 
beam  is  proportional  to  the  degree  to  which  optical  phase  conjugation  is  achieved.  The 
outgoing  field  after  reflection  from  DM2  is 

Uz(7)  =  U,(7)e^  .  (4.2) 

The  outgoing  field  after  reflection  from  DM1  is 

U0  (?)  =  Tz  [Uz  (r)]e*W  =  Tz  [u,  (7)e^r)]e^7)  .  (4.3) 

Again,  the  far-field  intensity  of  the  beam  is  proportional  to  the  degree  to  which  phase 
conjugation  is  achieved,  i.e., 

Ioc\\drUb{r)U0(rj\  .  (4.4) 

Substituting  Equation  (4.3)  into  Equation  (4.4)  and  applying  a  scale  factor  so  that  /  ranges 
from  0  to  1  yields 
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(4.5) 


J drU b  (¥)emTz  \u,  (r)eih<r)~\ 2 
\\drUb(r)U'b  {r)\\\dr,Ul{r')U]{r^ 

Now  the  goal  is  to  determine  phase  commands  that  maximize  /  when  applied  to  DM1  and 
DM2.  A  local  maximum  is  located  by  iteratively  solving  for  the  values  of  <f>i  and  fa  that 
maximize  I  (with  all  other  values  fixed).  Clearly,  the  numerator  of  Equation  (4.5)  is 
maximized  if  the  phase  of  the  laser  wavefront  is  the  conjugate  of  the  phase  of  the 
incoming  beacon  field.  This  condition  requires  the  integration  of  a  purely  real  number. 
With  all  other  values  fixed,  the  ideal  value  for  fa  is 

^  =  arg  jf/*  (r)Tx  [u,  (F>^  Jj  .  (4.6) 

An  expression  for  fa  is  found  with  the  aid  of  a  generalized  form  of  Rayleigh’s  theorem 
(i.e.,  the  energy  theorem,  Parseval’s  theorem,  or  Plancherel’s  theorem), 

\yaf{a)g(a)  =  \ypF{p)G'(P)  ,  (4.7) 

plus  a  basic  property  of  Fourier  transforms: 

f  (±x)  <=>  F*  (+^ )  .  (4.8) 

Applying  Equation  (4.7)  twice  and  invoking  Equation  (4.8)  to  simplify  results,  the 
formula  for  /  is 

|J  drU,  (r)eiMy  [ub 

1  ’  |J  drut  Fwi  (7|f  dfu,  (pju;  (rjf  •  <4'9) 

A  relationship  for  manipulating  these  expressions  is  formulated  using  Equation  (4.8)  and 
the  fact  that  for  energy  projection  applications  the  inversion  of  coordinates  is 
inconsequential.  Thus 
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or 


{7;  p/]}' =r[tr] 


(4.10) 


{^[t/*]}' =r[u\  .  (4,ii) 

From  Equation  (4.9)  the  optimal  value  of  <p2  with  all  other  quantities  fixed  is 

<t>2  =  arg \u]  (f)T\ Ub  (r  ) e* (P) T }  .  (4.12) 

/max  L  L  J 

The  above  algorithm  derivation  was  originally  developed  by  Ellerbroek  [11,12], 

4.2  Branch  Points  and  Least  Squares  Reconstruction 

Wavefront  sensors  normally  provide  data  that  indirectly  describe  the  wavefront 
phase.  In  the  case  of  the  Shack-Hartmann  sensor,  quantities  proportional  to  wavefront 
slopes  are  obtained.  The  slopes  are  related  to  the  phases  by 

g  =  GO  ,  (4.13) 

where  G  is  a  matrix  describing  the  geometrical  configuration  of  wavefront  sensor 
positions  and  phase  determination  positions  [4:258].  In  a  closed-loop  implementation, 
where  G  includes  the  effect  of  a  deformable  mirror,  G  is  often  referred  to  as  the  influence 
or  poke  matrix.  The  least  squares  estimate  for  the  phase  is 

6  =  (GrG)"1Gg  .'  (4.14) 

The  propagation  of  a  coherent  monochromatic  optical  field  through  a  turbulent 
atmosphere  (where  the  field  is  subject  to  spatially  varying  phase  perturbations)  results  in 
constructive  and  destructive  interference  in  the  amplitude  (or  intensity)  of  the  field. 
Locations  in  the  field  where  the  amplitude  goes  to  zero  indicate  the  presence  of  a  branch 
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point.  In  practice,  branch  points  in  an  optical  field  are  located  by  summing  the  principal 
value  gradients  around  the  smallest  possible  closed  contour,  where  the  principal  value 
operator  PF[]  simply  produces  an  equivalent  phase  in  the  range  -k  to  n.  If  a  branch 
point  is  enclosed,  then  the  contour  integral  of  principal  value  phase  gradients  equals 
±2 n  .  Likewise,  if  a  branch  point  is  not  enclosed  then  the  contour  integral  equals  zero. 
The  least  squares  estimator  for  the  phase  does  not  correctly  reconstruct  the  phase  if 
branch  points  are  present  because  of  the  underlying  assumption,  made  by  the  least 
squares  estimator,  that  the  measured  slopes  are  indicative  of  the  gradient  of  a  scalar  phase 
function.  Actually,  the  function  describing  the  gradient  of  the  perturbed  phase  must  be 
treated  as  the  sum  of  the  gradient  of  a  scalar  potential  and  the  curl  of  a  vector  potential, 
i.e., 

g(r)  =  Vs(r)+VxH(r)  ,  (4.15) 

where  r  is  a  position  vector  of  x  and  y  components,  s(r)  is  a  scalar  potential,  and  H(r)  is  a 
vector  potential.  Since  g(r)  only  has  components  in  the  x  andy  directions,  H(r)  clearly 
has  non-zero  components  only  in  the  z  direction. 

In  going  from  a  discrete  space  to  a  continuous  formulation,  the  G  matrix  above 
can  be  equated  to  the  gradient  operator,  GT  can  be  equated  to  the  divergence  operator, 
and  G  G  can  be  equated  to  the  laplacian  operator.  This  formulation  allows  additional 
insight  into  the  operation  of  the  least  squares  reconstructor.  Equation  (4. 14)  describing 
the  least  squares  estimate  is  recast  in  a  continuous  space  formulation  as 

V2€>(r)  =  V-[Vs(r)+VxH(r)]  .  (4.16) 
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Recalling  the  vector  identity  stating  that  the  divergence  of  the  curl  of  a  vector  is  equal  to 
zero,  Equation  (4.16)  becomes 

V24>  (r)  =  V  •  [V  s  (r)]  ,  (4.17) 

which  indicates  that  the  least  squares  estimate  for  the  phase  has  completely  ignored  the 
contribution  from  the  vector  potential.  The  actual  phase  can  now  be  defined  as 

0(r)  =  0(r)+0hid(O  >  (4.18) 

where  0hid(r)  is  the  “hidden”  phase  unaccounted  for  in  the  least  squares  reconstruction 
[13]. 

The  components  of  Equation  (4.18)  are  plotted  in  Figures  Figure  4. 1  through 
Figure  4.6  for  the  phase  of  one  realization  of  a  scintillated  test  field  with  a  theoretical  log- 
amplitude  variance  of  0.3.  Figure  4. land  Figure  4.2  show  a  bitmap  and  a  three- 
dimensional  surface  of  the  wrapped  phase,  respectively.  Figure  4.3  and  Figure  4.4  show 
the  least  squares  reconstruction  of  the  phase;  this  continuous  surface  is  a  form  that  is 
realizable  by  a  continuous  facesheet  deformable  mirror.  The  hidden  phase  (containing 
the  branch  points)  is  displayed  in  Figures  Figure  4.5  and  Figure  4.6  and  contains  abrupt 
changes  that  are  difficult  for  a  continuous  facesheet  deformable  mirror  to  realize  in 
practice. 
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Figure  4.1.  Wrapped  Phase  Map  for  One  Scintillated  Test  Field  with  a  Log- 
Amplitude  Variance  of  0.3  (Bitmap) 
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Figure  4.2.  Wrapped  Phase  Map  for  One  Scintillated  Test  Field  with  a  Log- 
Amplitude  Variance  of  0.3  (Surface) 
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Figure  4.3.  Least  Squares  Reconstruction  of  a  Wrapped  Phase  Map  for  One 
Scintillated  Test  Field  with  a  Log-Amplitude  Variance  of  0.3  (Bitmap) 


Figure  4.4.  Least  Squares  Reconstruction  of  a  Wrapped  Phase  Map  for  One 
Scintillated  Test  Field  with  a  Log-Amplitude  Variance  of  0.3  (Surface) 
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Pixels 

Figure  4.5.  Hidden  Phase  Contribution  of  a  Wrapped  Phase  Map  for  One 
Scintillated  Test  Field  with  a  Log-Amplitude  Variance  of  0.3  (Bitmap) 


Figure  4.6.  Hidden  Phase  Contribution  of  a  Wrapped  Phase  Map  for  One 
Scintillated  Test  Field  with  a  Log-Amplitude  Variance  of  0.3  (Surface) 
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5.  Results  and  Analysis 


Scintillated  beacon  fields  are  created  using  FFT-based  phase  screen  generation 
and  Fresnel  propagation  to  evaluate  the  performance  of  the  SGPA  algorithm.  All 
simulations  are  accomplished  using  MATLAB  version  5.3.  The  test  set  is  based  on  an 
ABL-type  scenario  and  North  Oscura  Peak  (NOP)  testing  in  support  of  ABL.  A  total  of 
32  fields  (and  the  corresponding  10  screen  atmospheric  models)  are  created  for  Rytov 
parameter  values  of  0.1,  0.2,. . .,  1.0.  The  outer  scale  ( L0 )  is  set  to  °°  for  all  simulations. 
Intensity  plots  of  one  scintillated  field  realization  for  each  Rytov  parameter  value  are 
provided  as  a  visual  reference  in  Appendix  A. 

For  comparison  the  phase  only  Strehl,  least  squares  phase  only  Strehl,  least 
squares  2-DM  Strehl,  and  uncompensated  Strehl  are  plotted  in  addition  to  the  2-DM 
Strehl  versus  the  Rytov  parameter.  Furthermore,  all  five  Strehls  are  computed  both 
inside  the  telescope  and  beyond  the  telescope  at  the  “top”  of  the  atmospheric  model  for 
the  Strehl  versus  Rytov  parameter  analysis.  The  “outer”  Strehl  is  computed  in  four 
different  versions  corresponding  to  increasing  levels  of  realism. 

5.1  Strehl  Ratio  vs.  the  Gaussian  Laser  Beam  Profile 

As  stated  previously,  the  input  aperture  of  the  telescope  is  set  to  31  pixels  in 
radius  for  all  simulations.  The  optimal  shape  of  the  Gaussian  amplitude  profile  is 
analyzed  by  plotting  Strehl  ratio  versus  location  of  the  1/e  point  in  the  Gaussian 
amplitude  field.  It  is  assumed  that  the  appropriate  optics  would  be  available  in  an  actual 
hardware  implementation  such  that  the  location  and  size  of  the  beam  waist  could  be 
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manipulated;  therefore,  the  shape  of  the  beam  is  optimized  in  the  plane  of  the  input 
aperture.  This  method  eliminates  any  dependence  of  the  optimized  value  on  the 
conjugate  range  of  DM2.  Also,  to  further  limit  the  influence  of  other  parameters  on  the 
results  of  this  simulation,  the  total  energy  in  the  Gaussian  beam  is  kept  constant  for  all 
beam  shapes.  This  is  easily  accomplished  by  examining  the  standard  equation  describing 
the  amplitude  of  a  Gaussian  beam  [14:70]: 


E(x,y,z ) 

wo 

—  _ avn 

f-T 

E0 

— - CAp 

w 

iw)_ 

where  w0  is  the  beam  waist  and  w  is  the  radial  location  of  the  1/  e  point  in  the  amplitude 
field  (which  is  a  function  of  distance  z  in  an  actual  Gaussian  beam).  For  the  purpose  of 
the  modeling  at  hand,  a  “default”  beam  waist  (w0)  is  arbitrarily  chosen  to  be  3 1  pixels  in 
radius.  Actually,  the  beam  waist  is  the  closest  the  1/e  point  ever  gets  to  the  axis  of 
propagation,  but  for  the  simulations  the  1/e  point  was  allowed  to  range  from  5  to  60 
pixels,  thereby  covering  nearly  all  of  the  128  x  128  beacon  field.  Setting  the  waist  to  a 
constant  simply  maintains  a  constant  energy  in  the  beam.  The  appropriate  scaling  of  the 
peak  amplitude  by  Equation  (5.1)  accomplishes  this  task  mathematically,  maintaining  a 
constant  value  of  the  integrated  intensity. 

Note  from  Figures  5.1  to  5.3  that  the  variation  of  Strehl  with  location  of  the  1/e 
point  decreases  with  increasing  Rytov  parameter.  From  this  perspective  stronger 
turbulence  is  advantageous  because  system  performance  is  less  dependent  on  the  degree 
and  accuracy  to  which  the  Gaussian  beam  profile  can  be  controlled.  The  peak  value  of 
the  curves  tends  to  increase  with  strength  of  turbulence  to  some  extent,  which  doesn’t 
correctly  describe  the  actual  performance  of  the  system  as  scintillation  increases.  The 
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next  step  toward  a  more  realistic  analysis  would  be  to  penalize  the  Strehl  for  laser  energy 
under  the  Gaussian  curve  not  accounted  for  by  the  finite  propagation  matrix,  which 
becomes  significant  as  the  waist  approaches  the  edge  of  the  matrix. 


Figure  5.1.  2-DM  Strehl  Ratio  (Inside  Telescope)  vs.  Gaussian  Laser  Beam  Waist 
Size  (in  Pixels)  for  Log-Amplitude  Variance  Values  of  0.1  -  0.5 


Figure  5.2.  2-DM  Strehl  Ratio  (Inside  Telescope)  vs.  Gaussian  Laser  Beam  Waist 
Size  (in  Pixels)  for  Log-Amplitude  Variance  Values  of  0.6  - 1.0 
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Figure  5.3.  Optimal  Gaussian  Laser  Beam  Waist  Size  in  the  Plane  of  the  Input 
Aperture  vs.  Log-Amplitude  Variance 

5.2  Strehl  Ratio  vs.  the  Conjugate  Range  of  the  2nd  DM 

Using  a  beam  waist  of  28  pixels  and  a  step-size  of  2  km,  the  mean  Strehl  (over  32 
realizations)  is  plotted  versus  the  conjugate  range  of  the  second  deformable  mirror  in 
Figure  5.4  and  Figure  5.5.  The  conjugate  range  corresponding  to  maximum  Strehl 
steadily  decreases  with  increasing  turbulence  strength  (Rytov  parameter),  which  is  clearly 
shown  in  Figure  5.6.  This  result  seems  non-intuitive  at  first  glance;  one  would  think  that 
the  opposite  might  be  true  because  a  greater  propagation  distance  provides  an  increased 
ability  of  phase  modulation  at  DM2  to  modulate  the  amplitude  profile  at  DM1 .  One 
possible  explanation  is  that  as  the  turbulence  strength  increases,  the  locations  of  regions 
of  significant  amplitude  in  the  sensed  beacon  field  vary  to  a  greater  extent.  This  effect 
requires  a  more  severe  modulation  of  the  phase  at  DM2,  which  in  turn  increases  the 
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energy  loss  due  to  the  finite  telescope  aperture.  The  optimal  range  must  therefore 
decrease  with  increasing  turbulence  in  order  to  overcome  the  energy  loss.  The  impact  of 
energy  loss  due  to  the  scattering  caused  by  DM2  was  first  discussed  by  Roggemann  [15]. 


Amplitude  Variance  Values  of  0.1  through  0.5 
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Optimal  Conjugate  Range  of  2nd  DM  (km)  In  2DM  Strehl  Ratio  (,nside  Telescope) 


Amplitude  Variance  Values  of  0.6  through  1.0 


Figure  5.6.  Optimal  Conjugate  Range  of  the  2nd  DM  vs.  Log-Amplitude  Variance 
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5.3  Number  of  Branch  Points  vs.  the  Conjugate  Range  of  the  2nd  DM 

The  conjugate  range  of  DM2  corresponding  to  the  maximum  number  of  branch 
points  in  DM1  steadily  decreases  with  increasing  turbulence  strength  but  saturates  near  5 
km  for  Rytov  values  of  0.5  and  greater;  see  Figure  5.7  and  Figure  5.8.  The  point  of 
saturation  resembles  that  of  the  actual  Rytov  versus  theoretical  Rytov  shown  in  Figure 
3.8.  Initial  inspection  of  Figure  5.9  and  Figure  5.10  indicates  a  much  larger  number  of 
branch  points  in  DM2  versus  DM  1 .  The  branch  points  in  DM2,  however,  are  counted 
over  the  entire  propagation  matrix  whereas  DM1  only  contains  branch  points  in  a  circular 
area  31  pixels  in  diameter;  recall  that  DM1  is  conjugate  to  the  collecting  aperture  of  the 
telescope  (which  acts  as  a  mask  on  DM1  during  iterations  of  the  SGPA).  The  sudden 
jumps  in  the  number  of  branch  points  in  DM2  are  due  to  the  fact  that  the  matrix 
dimension  is  allowed  to  jump  from  128  x  128  pixels  to  256  x  256  pixels  when  required  to 
avoid  propagator  aliasing.  The  overall  trend  in  the  number  of  branch  points  in  DM2  is 
decreasing  with  increasing  conjugate  range  of  DM2. 
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Figure  5.7.  Mean  Number  of  Branch  Points  in  DM1  vs.  the  Conjugate  Range  of  the 
2nd  DM  for  Log-Amplitude  Variance  Values  of  0.1  through  0.5 


Figure  5.8.  Mean  Number  of  Branch  Points  in  DM1  vs.  the  Conjugate  Range  of  the 
2nd  DM  for  Log-Amplitude  Variance  Values  of  0.6  through  1.0 
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Figure  5.9.  Mean  Number  of  Branch  Points  in  DM2  vs.  the  Conjugate  Range  of  the 
2nd  DM  for  Log-Amplitude  Variance  Values  of  0.1  through  0.5 


Figure  5.10.  Mean  Number  of  Branch  Points  in  DM2  vs.  the  Conjugate  Range  of 
the  2nd  DM  for  Log-Amplitude  Variance  Values  of  0.6  through  1.0 
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5.4  Number  of  Iterations  to  Convergence  vs.  the  Conjugate  Range  of  the  2nd  DM 


The  mean  number  of  iterations  to  convergence  of  the  SGPA  (to  within  0.0005)  is 
plotted  versus  the  conjugate  range  of  the  second  DM  in  Figure  5.1 1  and  Figure  5.12.  The 
peak  number  of  iterations  per  Rytov  value  varies  between  approximately  60  and  80  and 
the  location  of  the  peak  decreases  with  increasing  turbulence  strength.  The  general 
shapes  of  the  curves  are  very  similar  to  the  shapes  describing  the  number  of  branch  points 
in  DM1 ;  see  Figure  5.7  and  Figure  5.8.  This  would  seem  to  indicate  that  a  larger  number 
of  SGPA  iterations  results  in  a  larger  number  of  branch  points  in  DM1 . 


Figure  5.11.  Mean  Number  of  Iterations  to  Convergence  vs.  the  Conjugate  Range  of 
the  2nd  DM  for  Log-Amplitude  Variance  Values  of  0.1  through  0.5 
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Figure  5.12.  Mean  Number  of  Iterations  to  Convergence  vs.  Conjugate  Range  of  the 
2nd  DM  for  Log-Amplitude  Variance  Values  of  0.6  through  1.0 

5.5  Strehl  Ratio  vs.  the  Rytov  Parameter 

Using  a  finite  input  aperture  of  3 1  pixels  and  a  beam  waist  size  (in  the  plane  of 
the  input  aperture)  of  28  pixels,  the  mean  Strehl  ratio  is  determined  over  an  ensemble  of 
32  realizations  for  each  value  of  the  Rytov  parameter.  Both  DMs  are  assumed  to  be 
“infinite”  in  size  for  this  simulation.  The  conjugate  range  of  the  2nd  DM  is  set  to  50  km 
except  for  the  case  of  cx2  =  1 .0  in  which  a  conjugate  range  of  48  km  is  used  due  to 
limited  computing  resources.  To  avoid  propagator  aliasing  ranges  greater  than  48  km 
require  a  matrix  size  of  512  x  512  or  greater  for  <3X2  =  1.0,  which  increases  the  required 
computing  time  significantly. 

The  Strehl  ratios  computed  inside  the  telescope  are  plotted  versus  the  Rytov 
parameter  in  Figure  5.13,  including  the  2-DM  Strehl,  2-DM  least  squares  Strehl,  1-DM 
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(i.e.,  phase  only)  Strehl,  1-DM  least  squares  Strehl,  and  the  uncompensated  Strehl.  For 
the  Strehl  ratios  computed  inside  the  telescope,  the  2-DM  Strehl  remained  above  0.9  for 
all  values  of  log-amplitude  variance.  The  performance  of  the  1-DM  system  gradually 
decreases  from  a  Strehl  of  0.75  to  slightly  less  than  0.6  for  Rytov  values  of  0.8  through 
1.0. 


Log-Amplitude  Variance  (a^) 


Figure  5.13.  Mean  Strehl  Ratio  Inside  Telescope  vs.  Log-Amplitude  Variance 


In  order  to  provide  a  more  realistic  evaluation  of  the  performance  of  the  2  DM 
AO  system,  the  Strehl  ratio  is  also  calculated  outside  the  telescope  by  propagating  the 
pre-compensated  outgoing  laser  beam  back  through  the  original  set  of  phase  screens. 
This  “back-propagated”  Strehl  is  calculated  in  several  ways  (to  provide  increasing  levels 
of  realism)  and  is  plotted  in  Figures  5.14  to  5.17. 
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First,  the  back-propagated  Strehl  is  calculated  using  the  original  laser  profile  as 
the  “ideal”  field.  Second,  the  laser  field  is  propagated  to  the  top  of  the  atmosphere  for 
use  as  the  “ideal”  field.  Third,  the  laser  is  propagated  to  the  telescope  aperture,  through 
the  aperture,  then  to  the  top  of  the  atmosphere.  Finally,  the  most  realistic  approach  is 
identical  to  the  previous  one  with  one  exception;  the  Strehl  ratio  is  penalized  for  energy 
lost  due  to  the  finite  aperture.  As  expected,  this  last  approach  results  in  the  lowest  Strehl 
ratios.  All  2-DM  Strehl  ratios  calculated  beyond  the  telescope,  however,  follow  a  similar 
trend  and  decrease  monotonically  with  increasing  turbulence  strength.  The  system 
performance  predicted  by  this  assessment  is  significantly  worse  than  that  predicted  by  the 
in-telescope  Strehls.  However,  even  the  most  critical  analysis  indicates  a  performance 
increase  of  approximately  0.1  Strehl  over  the  1-DM  system  for  the  worst  turbulence 
scenario. 


Ul  Not  Propagated 


Figure  5.14.  Mean  Strehl  Ratio  Outside  the  Telescope  vs.  Log-Amplitude  Variance 
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Mean  Strehl  (Outside  Telescope)  ^  Mean  Strehl  (Outside  Telescope) 


Ui  Propagated 


'.15.  Mean  Strehl  Ratio  Outside  Telescope  vs.  Log-Amplitude  Variance 


UI  Propagated  Through  Aperture 


Log-Amplitude  Variance  (o^) 


Figure  5.16.  Mean  Strehl  Ratio  Outside  Telescope  vs.  Log-Amplitude  Variance 
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Ul  Propagated  Though  Aperture  (with  Penalty  for  Energy  Loss) 


Figure  5.17.  Mean  Strehl  Ratio  Outside  Telescope  vs.  Log-Amplitude  Variance 


Figures  5. 1 8  to  5.22  plot  the  different  Strehl  ratios  by  type  of  AO  system  using 
both  the  in-telescope  and  beyond-telescope  approaches.  For  all  systems  the  first  and 
second  methods  for  calculating  the  Strehl  ratio  beyond  the  telescope  yield  nearly  identical 
results.  This  result  is  not  surprising  given  the  fact  that  the  Fourier  transform  of  a 
Gaussian  is  another  Gaussian.  Furthermore,  the  divergence  angle  characterizing  the 
Gaussian  beam  is  not  incorporated  into  this  analysis. 
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Mean  2  DM  LS  Strehl  Ratio 


Figure  5.19.  Mean  Least  Squares  2-DM  Strehl  Ratio  vs.  Log-Amplitude  Variance 
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Mean  1  DM  LS  Strehi  Ratio 


Figure  5.21.  Mean  Least  Squares  1-DM  Strehi  Ratio  vs.  Log-Amplitude  Variance 
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Figure  5.22.  Mean  Uncompensated  Strehl  Ratio  vs.  Log-Amplitude  Variance 


5.6  Number  of  Branch  Points  vs.  the  Rytov  Parameter 

The  numbers  of  branch  points  in  each  of  the  deformable  mirror  commands  are 
plotted  versus  the  log-amplitude  variance  in  Figure  5.23  and  Figure  5.24.  The  number  of 
branch  points  in  the  beacon  field  or,  conversely,  in  the  1-DM  system  mirror  commands 
are  plotted  versus  Rytov  parameter  in  Figure  5.25.  The  number  of  branch  points  in  the 
DM1  commands  for  the  2-DM  system  is  nearly  600  branch  points  for  the  weakest 
turbulence,  but  decreases  to  about  240  branch  points  for  Rytov  values  of  0.9  and  1 .0.  The 
numbers  of  branch  points  in  the  DM2  commands  are  orders  of  magnitude  greater  for  log- 
amplitude  variance  values  above  0.5.  Once  again,  this  is  due  to  counting  branch  points 
over  a  larger  grid;  a  256  x  256  pixel  propagation  matrix  is  required  for  Rytov  values 
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greater  than  or  equal  to  0.6  and  propagation  distances  greater  than  or  equal  to  50  km. 
Furthermore,  the  branch  points  in  the  DM1  commands  are  effectively  limited  to  the  area 
within  a  circular  aperture  of  radius  equal  to  31  pixels  due  to  the  fact  that  DM1  is 
conjugate  to  the  collecting  aperture  of  the  telescope.  The  large  numbers  of  branch  points 
explains  the  extremely  poor  performance  (only  slightly  better  than  uncompensated)  of  the 
least  squares  2-DM  system,  as  seen  in  Figure  5.19.  The  hidden  phase  for  these  numbers 
of  branch  points  clearly  contains  a  significant  portion  of  the  information  for  the  DM 
controls.  The  relatively  small  number  of  branch  points  in  the  1-DM  system  mirror 
commands  (i.e.,  beacon  field)  is  one  reason  that  phase-only  correction  is  able  to  maintain 
significant  performance  increases  over  an  uncompensated  system  even  at  the  worst 
turbulence  strengths. 


Figure  5.23.  Mean  Number  of  Branch  Points  in  DM1  vs.  Log-Amplitude  Variance 


5-19 


0  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9  1 


Log-Amplitude  Variance  (cr^) 

Figure  5.24.  Mean  Number  of  Branch  Points  in  DM2  vs.  Log-Amplitude  Variance 
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Figure  5.25.  Mean  Number  of  Branch  Points  in  Beacon  Field  (i.e.,  1-DM  System 
Mirror  Commands)  vs.  Log-Amplitude  Variance 
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5. 7  Number  of  Iterations  to  Convergence  vs.  the  Rytov  Parameter 

The  number  of  iterations  of  the  SGPA  algorithm  required  to  converge  to  within 
0.0005  Strehl  decreases  rapidly  with  increasing  Rytov  parameter  up  to  a  Rytov  value  of 
0.5;  see  Figure  5.26.  The  saturation  of  the  required  iterations  near  27  occurs  at  about  az2 
=  0.5,  which  roughly  corresponds  to  the  value  of  theoretical  Rytov  past  which  a 
saturation  in  measured  Rytov  (see  Figure  3.8)  occurs.  The  decrease  in  required  iterations 
may  be  due  to  the  fact  that  as  scintillation  becomes  worse  there  is  less  integrated 
amplitude  in  the  beacon  field  (to  which  the  SGPA  must  reshape  the  Gaussian  laser  beam 
against). 


Figure  5.26.  Mean  Number  of  Iterations  to  SGPA  Convergence  vs.  Log-Amplitude 

Variance 


5-21 


5.8  Strehl  Ratio  vs.  the  Radii  of  the  Deformable  Mirrors 


To  model  the  impact  of  energy  loss  (due  to  finite  DM  size)  on  the  Strehl  ratio,  the 
radii  of  both  DMs  are  varied  simultaneously.  This  is  accomplished  in  two  ways:  1)  by 
masking  the  field  (not  just  the  phase),  which  is  analogous  to  forming  a  finite  mirror  for 
which  the  entire  surface  is  deformable  and  2)  by  masking  the  phase,  which  is  analogous 
to  forming  an  infinite  flat  surface  mirror  for  which  only  a  finite  center  portion  is 
deformable.  A  real  system  would  probably  mask  the  phase  (i.e.,  use  a  large  DM  but  only 
control  a  central  section).  However,  the  results  presented  here  indicate  only  a  very  minor 
difference  in  Strehl  for  the  two  approaches.  Both  approaches  are  used  to  plot  the  2-DM 
Strehl,  the  least  squares  2-DM  Strehl,  the  1-DM  Strehl,  and  the  least  squares  1-DM  Strehl 
in  Figure  5.27  for  a  Rytov  value  of  0.1.  The  results  are  so  similar,  in  fact,  that  only  the 
first  (most  conservative)  approach  is  used  for  the  remaining  simulations  in  this  section. 


Figure  5.27.  Strehl  Ratio  (Inside  Telescope)  vs.  Radii  of  the  Deformable  Mirrors  for 

a  Log-Amplitude  Variance  of  0.1 
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The  SGPA  is  run  once  for  each  beacon  field  realization  per  Rytov  value,  and  thus 
the  subsequent  Strehl  ratios  (for  a  given  beacon  field  realization)  are  all  calculated  with 
one  set  of  phase  commands  with  varying  mirror  sizes.  The  results  are  plotted  in  Figures 
5.28  to  5.37.  As  expected,  the  Strehl  of  the  1-DM  system  reaches  an  absolute  maximum 
as  the  radius  of  the  single  DM  equals  the  radius  of  the  input  aperture.  This  outcome  is 
due  to  the  fact  that  the  1-DM  system  can  only  impact  the  outgoing  field  in  the  plane  of 
the  telescope’s  collecting  aperture.  The  2-DM  system,  however,  achieves  a  significant 
increase  in  Strehl  for  DM  radii  that  increase  beyond  the  size  of  the  input  aperture.  This 
result  also  makes  sense  due  to  the  fact  that  phase  modulations  at  any  point  in  the  outgoing 
field  in  the  plane  of  DM2  can  affect  the  amplitude  of  the  field  at  (theoretically)  any  other 
point  in  the  field  after  propagation  to  DM1 . 


Figure  5.28.  Strehl  Ratio  (Inside  Telescope)  vs.  Radii  of  the  Deformable  Mirrors  for 

a  Log-Amplitude  Variance  of  0.1 


5-23 


Figure  5.29.  Strehl  Ratio  (Inside  Telescope)  vs.  Radii  of  the  Deformable  Mirrors  for 

a  Log-Amplitude  Variance  of  0.2 


Figure  5.30.  Strehl  Ratio  (Inside  Telescope)  vs.  Radii  of  the  Deformable  Mirrors  for 

a  Log-Amplitude  Variance  of  0.3 


5-24 


Figure  5.31.  Strehl  Ratio  (Inside  Telescope)  vs.  Radii  of  the  Deformable  Mirrors  for 

a  Log- Amplitude  Variance  of  0.4 


Radii  of  Deformable  Mirrors  (Pixels) 


Figure  5.32.  Strehl  Ratio  (Inside  Telescope)  vs.  Radii  of  the  Deformable  Mirrors  for 

a  Log-Amplitude  Variance  of  0.5 


5-25 


o2  =  0.6 


Figure  5.33.  Strehl  Ratio  (Inside  Telescope)  vs.  Radii  of  the  Deformable  Mirrors  for 

a  Log-Amplitude  Variance  of  0.6 


Figure  5.34.  Strehl  Ratio  (Inside  Telescope)  vs.  Radii  of  the  Deformable  Mirrors  for 

a  Log-Amplitude  Variance  of  0.7 
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o2  =  0.8 


Figure  5.35.  Strehl  Ratio  (Inside  Telescope)  vs.  Radii  of  the  Deformable  Mirrors  for 

a  Log-Amplitude  Variance  of  0.8 


Figure  5.36.  Strehl  Ratio  (Inside  Telescope)  vs.  Radii  of  the  Deformable  Mirrors  for 

a  Log-Amplitude  Variance  of  0.9 
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Radii  of  Deformable  Mirrors  (Pixels) 


Figure  5.37.  Strehl  Ratio  (Inside  Telescope)  vs.  Radii  of  the  Deformable  Mirrors  for 

a  Log-Amplitude  Variance  of  1.0 


5.9  Number  of  Branch  Points  vs.  the  Radii  of  the  Deformable  Mirrors 

The  numbers  of  branch  points  in  each  of  the  DMs  is  plotted  versus  log-amplitude 
variance  in  Figures  5.38  to  5.47.  As  with  the  phase-only  correction  system,  the  numbers 
of  branch  points  in  the  controls  for  DM1  reach  an  absolute  maximum  at  DM  radii  equal 
to  and  greater  than  that  of  the  input  aperture.  This  result  is  due  to  the  fact  that  the  first 
deformable  mirror  can  only  effect  the  phase  of  the  outgoing  laser  in  the  area  limited  by 
the  collecting  aperture  of  the  telescope.  This  fact  is  expressed  mathematically  in  the 
SGPA  phase  iterative  update  equation  for  DM1  commands  in  Equation  (4.6).  The 
numbers  of  branch  points  in  the  commands  for  DM2  undergo  a  sharp  increase  for  radii 
near  45  pixels  for  all  levels  of  turbulence  strength. 
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Figure  5.38.  Number  of  Branch  Points  in  2-DM  Commands  vs.  Radii  of  DMs  for  a 

Log-Amplitude  Variance  of  0.1 


Figure  5.39.  Number  of  Branch  Points  in  2-DM  Commands  vs.  Radii  of  DMs  for  a 

Log-Amplitude  Variance  of  0.2 
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o2  =  0.3 


Figure  5.40.  Number  of  Branch  Points  in  2-DM  Commands  vs.  Radii  of  DMs  for  a 

Log-Amplitude  Variance  of  0.3 


a2  =  0.4 


Figure  5.41.  Number  of  Branch  Points  in  2-DM  Commands  vs.  Radii  of  DMs  for  a 

Log-Amplitude  Variance  of  0.4 
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o2  =  0.5 


Figure  5.42.  Number  of  Branch  Points  in  2-DM  Commands  vs.  Radii  of  DMs  for  a 

Log- Amplitude  Variance  of  0.5 


ct2  =  0.6 


Figure  5.43.  Number  of  Branch  Points  in  2-DM  Commands  vs.  Radii  of  DMs  for  a 

Log-Amplitude  Variance  of  0.6 
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Figure  5.44.  Number  of  Branch  Points  in  2-DM  Commands  vs.  Radii  of  DMs  for  a 

Log-Amplitude  Variance  of  0.7 
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Figure  5.45.  Number  of  Branch  Points  in  2-DM  Commands  vs.  Radii  of  DMs  for  a 

Log-Amplitude  Variance  of  0.8 
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Radii  of  Deformable  Mirrors  (Pixels) 


Figure  5.46.  Number  of  Branch  Points  in  2-DM  Commands  vs.  Radii  of  DMs  for  a 

Log-Amplitude  Variance  of  0.9 


Figure  5.47.  Number  of  Branch  Points  in  2-DM  Commands  vs.  Radii  of  DMs  for  a 

Log-Amplitude  Variance  of  1.0 
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6.  Conclusions  and  Recommendations 


The  addition  of  a  second  deformable  mirror  to  the  two-DM  AO  transmission 
system  clearly  improves  theoretical  performance  (over  a  single  DM  system)  for 
delivering  energy  on  target  through  atmospheric  turbulence.  Strehl  ratios  calculated  by 
back-propagating  the  modulated  laser  beam  through  the  atmosphere,  however,  indicate  a 
less  significant  improvement  over  the  single  DM  system  (compared  to  Strehl  ratios 
calculated  inside  the  telescope). 

The  hidden  phase  contained  in  the  branch  points  is  critical  to  the  performance  of 
the  SGPA,  as  indicated  by  the  two-DM  least  squares  Strehl  ratio  results.  The  DM 
commands  generated  by  the  SGPA  are  corrupted  by  a  large  number  of  branch  points,  and 
this  problem  has  at  least  two  possible  solutions.  One  solution  is  the  use  of  a  phase 
reconstruction  scheme  other  than  least  squares,  and  another  solution  is  the  use  of  a  new 
two-DM  algorithm  that  imposes  constraints  on  the  number  of  branch  points  in  the  DM 
commands.  An  effective  real-world  solution  will  probably  combine  both  approaches,  and 
work  is  underway  toward  this  end. 

Allowing  the  conjugate  range  of  DM2  to  be  finite  minimizes  the  energy  lost  at 
DM1  due  to  scattering  caused  by  the  phase  modulation  at  DM2.  For  this  reason  the  two- 
DM  system  with  the  second  DM  conjugate  to  a  finite  range  should  outperform  the  system 
with  the  second  DM  conjugate  to  the  far-field  (infinity). 

Recommendations  for  further  research  are  provided  below. 
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First,  track  computing  time  per  SGPA  iteration.  Such  tracking  could  be  done  for 
all  of  the  analyses  accomplished  in  this  thesis,  with  the  final  result  being  an  average 
number  of  seconds  per  iteration. 

Second,  perform  the  Strehl  versus  Gaussian  beam  waist  size  optimization  with  the 
inclusion  of  a  penalty  to  the  Strehl  ratio  calculation  for  energy  in  the  Gaussian  laser  beam 
lost  (i.e.,  unaccounted  for)  outside  the  finite  propagation  matrix. 

Third,  execute  the  Strehl  versus  DM  radii  analysis  while  keeping  track  of  the 
amount  of  energy  lost  outside  the  finite  DMs.  The  amount  of  lost  energy  should  increase 
with  decreasing  DM  size  and  with  increasing  turbulence  strength.  Similarly,  the  Strehl 
versus  conjugate  range  of  DM2  analysis  could  be  accomplished  while  keeping  track  of 
the  amount  of  energy  lost  outside  of  DM1 .  The  energy  lost  will  most  likely  increase  with 
increasing  conjugate  range  of  DM2  and  increasing  turbulence  strength. 

Finally,  complete  the  number  of  branch  points  in  DM2  versus  both  Rytov  value 
and  the  conjugate  range  of  DM2  analyses  while  constraining  the  area  within  which 
branch  points  are  counted  to  the  size  of  the  finite  collecting  aperture  of  the  telescope. 

This  extension  would  provide  a  basis  for  fair  comparison  of  the  number  of  branch  points 
in  DM1  with  the  number  of  branch  points  in  DM2,  due  to  the  fact  that  the  numbers  of 
branch  points  in  DM1  are  inherently  counted  only  within  an  area  equal  to  the  size  of  the 
collecting  aperture,  as  DM1  is  conjugate  to  that  plane.  Another  approach  would  be  to 
normalize  the  number  of  branch  points  to  some  unit  area.  The  locations  of  branch  points 
could  also  be  investigated.  It  has  been  reported  in  the  literature  that  branch  points  tend  to 
occur  more  often  in  regions  of  lower  intensity  [12]. 
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In  summary,  the  more  critical  Strehl  ratio  analysis  performed  by  back-propagating 
the  pre-compensated  laser  beam  through  the  atmospheric  models  indicates  that  a  two-DM 
AO  system  provides  increased  performance  over  that  of  a  one-DM  system.  Accounting 
for  the  hidden  phase  contained  in  the  branch  points  corrupting  the  mirror  commands 
(produced  by  the  SGPA  algorithm),  however,  is  critical  to  successful  implementation  of 
the  SGPA. 
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Appendix  A.  Fields  Used  for  Testing  SGPA  Algorithm 


Pixels 

Figure  A.l.  Intensity  Plot  of  One  Realization  for  a  Log-Amplitude  Variance  of  0.1 


Pixels 

Figure  A.2.  Intensity  Plot  of  One  Realization  for  a  Log-Amplitude  Variance  of  0.2 
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Pixels 

Figure  A.3.  Intensity  Plot  of  One  Realization  for  a  Log- Amplitude  Variance  of  0.3 
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Figure  A.4.  Intensity  Plot  of  One  Realization  for  a  Log-Amplitude  Variance  of  0.4 
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Figure  A.5.  Intensity  Plot  of  One  Realization  for  a  Log-Amplitude  Variance  of  0.5 
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Figure  A.6.  Intensity  Plot  of  One  Realization  for  a  Log-Amplitude  Variance  of  0.6 
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Figure  A.7.  Intensity  Plot  of  One  Realization  for  a  Log-Amplitude  Variance  of  0.7 
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Figure  A.8.  Intensity  Plot  of  One  Realization  for  a  Log-Amplitude  Variance  of  0.8 
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Figure  A.9.  Intensity  Plot  of  One  Realization  for  a  Log- Amplitude  Variance  of  0.9 
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Figure  A.10.  Intensity  Plot  of  One  Realization  for  a  Log-Amplitude  Variance  of  1.0 
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